function [betas, irtc, ptc, se] = calcbetas(rtc, tcd, dim, tol)
% calcbetas  - perform GLM regression
%
% FORMAT:       [betas, irtc, ptc, se] = calcbetas(X, tcd [, tdim [, tol]])
%
% Input fields:
%
%       X           design matrix (TxR)
%       tcd         time course data (N-dim)
%       tdim        temporal dimension (default: 1)
%       tol         tolerance value to set iXX to 0 (default: 4 * eps)
%
% Output fields:
%
%       betas       betas, in time course data dimension
%       irtc        inverse design matrix
%       ptc         predicted tc
%       se          standard error

% Version:  v0.8a
% Build:    9102316
% Date:     Oct-23 2009, 4:44 PM CEST
% Author:   Jochen Weber, SCAN Unit, Columbia University, NYC, NY, USA
% URL/Info: http://wiki.brainvoyager.com/BVQXtools

% argument check
if nargin < 2 || ...
   ~isa(rtc, 'double') || ...
    ndims(rtc) ~= 2 || ...
    isempty(rtc) || ...
    any(isinf(rtc(:)) | isnan(rtc(:))) || ...
   ~isnumeric(tcd) || ...
   ~any(size(tcd) == size(rtc, 1)) || ...
    isempty(tcd)
    error( ...
        'BVQXtools:BadArgument', ...
        'Bad design matrix or time course data supplied.' ...
    );
end

% further arguments
if nargin < 3 || ...
   ~isa(dim, 'double') || ...
    isempty(dim)
    rdim = findfirst(size(tcd) == size(rtc, 1));
elseif numel(dim) ~= 1 || ...
    isinf(dim) || ...
    isnan(dim) || ...
    fix(dim) ~= dim || ...
    dim < 1 || ...
    dim > 5
    error( ...
        'BVQXtools:BadArgument', ...
        'Bad tdim argument supplied.' ...
    );
else
    rdim = real(dim);
end
if nargin < 4 || ...
   ~isa(tol, 'double') || ...
    isempty(tol)
    tol = 4 * eps;
elseif numel(tol) ~= 1 || ...
    isnan(tol) || ...
    abs(tol) > sqrt(eps)
    error( ...
        'BVQXtools:BadArgument', ...
        'Bad tol argument supplied.' ...
    );
else
    tol = abs(tol);
end

% build rtc' and pinv(rtc' * rtc)
try
    trtc = rtc';
    irtc = pinv(trtc * rtc);
    irtc(abs(irtc) < tol) = 0;
    prtc = irtc * trtc;
catch
    error( ...
        'BVQXtools:MATLABError', ...
        'Error building (pseudo) inverse of design matrix: %s.', ...
        lasterr ...
    );
end
    
% number of rows and betas
numrows = size(rtc, 1);
numbets = size(rtc, 2);

% make data double
if ~isa(tcd, 'double')
    if istransio(tcd)
        tcd = double(resolve(tcd));
    else
        tcd = double(tcd);
    end
end

% check dim
if rdim < 1 || ...
    rdim > length(size(tcd)) || ...
    size(tcd, rdim) ~= numrows
    error( ...
        'BVQXfile:InternalError', ...
        'Invalid data matrix vs. selected dimension.' ...
    );
end
if rdim > 1
    neword = [rdim, 1:length(size(tcd))];
    newodr = find(neword == rdim);
    neword(newodr(2)) = [];
    tcd = permute(tcd, neword);
    [newsrt, oldord] = sort(neword);
else
    oldord = [];
end

% reshaping data to comply
tcds = size(tcd);
numvox = prod(tcds(2:end));
tcd = reshape(tcd, [tcds(1), numvox]);
tcds(1) = [];
tcdrs = tcds;
if length(tcdrs) < 2
    tcdrs(2) = 1;
end

% perform calculus
betas = prtc * tcd;

% more output?
if nargout > 2
    ptc = rtc * betas;
    if nargout > 3
        se = squeeze(reshape(std(tcd - ptc, 0), tcdrs)) ...
             * sqrt((numrows - 1) / (numrows - numbets));
    end
    ptc = reshape(ptc, [numrows, tcds]);
    if ~isempty(oldord)
        ptc = permute(ptc, oldord);
    end
end

% reshape betas
betas = reshape(shiftdim(betas, 1), [tcds, numbets]);
